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i!B ■ ABSTRACT 

o , 

^ , We investigate the growth of spiral-arm substructure in vertically stratified, 

^ I self-gravitating, galactic gas disks, using local numerical MHD simulations. Our 

new models extend our previous two-dimensional studies (Kim & Ostriker 2002), 
I which showed that a magnetized spiral shock in a thin disk can undergo magneto- 

^ I Jeans instability (MJI), resulting in regularly-spaced interarm spur structures 

and massive gravitationally-bound fragments. Similar spur (or "feather") fea- 
tures have recently been seen in high-resolution observations of several galaxies, 
and massive bound gas condensations are likely the precursors of giant molec- 
ular cloud complexes (GMCs) and H II regions. Here, we consider two sets of 
numerical models: two-dimensional simulations that use a "thick-disk" gravita- 
tional kernel, and three-dimensional simulations with explicit vertical stratifica- 
tion. Both models adopt an isothermal equation of state with Cg = 7 km s~^. 
When disks are sufficiently magnetized and self-gravitating, the result in both 
sorts of models is the growth of spiral-arm substructure similar to that in our 
previous razor-thin models. Reduced self-gravity due to nonzero disk thickness 
increases the spur spacing to ~ 10 times the Jeans length at the arm peak, 
a factor ~ 3 — 5 times larger than in razor-thin models. Bound clouds that 
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form from spur fragmentation have masses ~ (1 — 3) x 10^ Mq each, a factor 
~ 3 — 8 times larger than in razor-thin models with the same gas surface density 
and stellar spiral arm strength. These condensation masses are comparable to 
results from other three-dimensional models without spiral structure, and simi- 
lar to the largest observed GMCs. The mass-to-fiux ratios and specific angular 
momenta of the bound condensations are lower than large-scale galactic values, 
as is true for observed GMCs. We find that unmagnetized or weakly magne- 
tized two-dimensional models are unstable to the "wiggle instability" previously 
identified by Wada & Koda (2004), and proposed as a potential spur- and clump- 
forming mechanism. However, our fully three-dimensional models do not show 
this effect. Non-steady motions and strong vertical shear prevent coherent vor- 
tical structures from forming, evidently suppressing the wiggle instability that 
appears in two-dimensional (isothermal) simulations. We also find no clear traces 
of Parker instability in the nonlinear spiral arm substructures that emerge (in 
self-gravitating models), although conceivably Parker modes may help seed the 
MJI at early stages since azimuthal wavelengths are similar. 

Subject headings: galaxies: ISM — instabilities — ISM: kinematics and dynamics 
— ISM: magnetic fields — MHD — stars: formation 

1. Introduction 

The interstellar medium (ISM) contains structure at many scales, and the growth of 
much of this structure is believed to represent the first step in initiating galactic star forma- 
tion. In particular, the molecular portion of the ISM is concentrated into dense clouds over 
a range of scales (Blitz 1993); each of these clouds also contains significant higher- density 
substructure. The largest of these cloud complexes (up to ~ 10^ Mq, including atomic 
envelopes) contain most of the mass in the overall distribution of giant molecular clouds 
(GMCs). In galaxies with prominent spiral structure, GMCs tend to be concentrated into 
the spiral arms. Observationally, the higher-than-average specific star formation rates (i.e. 
per unit gas mass) within spiral arms (Knapen et al. 1992, 1996) suggests an evolution- 
ary sequence: compression of diffuse gas (and/or collection of smaller clouds) within the 
arms prompts GMC formation, and star formation is subsequently triggered in the densest 
portions within GMCs. 

Prom a theoretical point of view, interaction of the ISM with spiral arms should naturally 
result in bound cloud formation when the mean density is high enough that self-gravitating 
instabilities are able to grow. The growth of self-gravitating structures in spiral arms has 



-3- 



been investigated by many authors, using both hnear-theory analyses (Balbus & Cowie 
1985; Balbus 1988; Elmegreen 1994; Kim & Ostriker 2002) and more recently, nonlinear 
hydrodynamic and magnetohydrodynamic (MHD) simulations that yield bound clouds with 
properties similar to observed GMCs (Kim & Ostriker 2002, hereafter Paper I). 

Interestingly, self-gravity acting on the gas within spiral arms can apparently also lead 

to growth of structures larger than GMCs. These structures take the shape of gaseous 
spurs that project from the main body of the arm into the interarm region. The first direct 
numerical simulations showing the development of these trailing spurs (Paper 1) occurred 
contemporaneously with the release of the now-famous Hubble Heritage image of M51 (Scov- 
ille & Rector 2001; Scoville et al. 2001), in which strikingly similar structures are evident. 
Furthermore, a recent Spitzer Legacy image of M51 displays a quasi-regular distribution of 
thin, traihng dust filaments throughout the interarm regions (Kennicutt 2004). The dust 
filaments seen in these images could well be the remnants of gaseous spurs initiated inside 
spiral arms, which have maintained their integrity deep into the interarm regions. A recent 
archival study of HST galaxy images (Lavigne, Vogcl, & Ostriker 2006, in preparation) has 
shown that arm substructures of the kind seen in M51 are in fact ubiquitous in the grand 
design spirals where the global pattern is clean enough to make identification of spiral spur 
substructures possible. When identified as a regular series of dust lanes extending out of 
a primary dust lane across the bright stellar arm and into the interarm region, such spur 
structures have historically been termed "feathers" (Lynds 1970). 

In this work, we extend the (thin disk) two-dimensional numerical models of Paper I into 
three dimensions, in order to study the effects of vertical stratification on the development of 
condensation instabilities in spiral arms. In our previous work, we argued that self-gravity, 
aided by magnetic fields, is the key to gaseous spur and dense cloud formation. However, 
it has recently been proposed that Kelvin- Helmholtz instabihties are also able to trigger 
spiral arm spur formation (Wada & Koda 2004). In addition, a theoretical proposal with a 
long history is that Parker instabilities in spiral arms may be important in prompting GMC 
formation (e.g. Parker 1966; Mouschovias, Shu, & Woodward 1974; Blitz & Shu 1980; see 
discussion and references in Kim, Ostriker, & Stone 2002, hereafter Paper 11). Our three- 
dimensional models allow us to explore both of these proposals using more realistic numerical 
set-ups than in previous work. 

We shall show that, for vertically stratified magnetized disks, qualitatively quite similar 
results to the conclusions of Paper 1 appear to hold. Quantitatively, we shall also show 
that the three-dimensional results are consistent with results of two-dimensional models 
when "thick-disk" dilution of gravity is incorporated. We also find that three-dimensional 
dynamics tends to suppress the growth of Kelvin- Helmholtz-driven spurs described by Wada 
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& Koda (2004) based on unmagnetized two-dimensional simulations, and that there are no 
clear signatures of Parker instability in our models. 

Star formation takes place within molecular clouds, and at least the most massive GMCs 
that dominate the distribution appear to live for significantly less than a galactic orbital time, 
based on their close association with spiral arms (Solomon, Sanders, & Rivolo 1985; Solomon 
& Rivolo 1989; Engargiola et al. 2003; Stark & Lee 2005b). Thus, the galactic star formation 
rate at a given radius can be written as a cloud formation rate Mdoud times a star 
formation efficiency per cloud lifetime, csf- To the extent that this star formation efficiency 
depends primarily on physical processes within a cloud, understanding the regulation of star 
formation can therefore be separated into a large-scale (galaxy) and a small-scale (GMC) 
problem. The current work primarily addresses the large-scale problem by investigating the 
mechanisms and timescales required for diffuse gas to be gathered into bound clouds in the 
case when strong spiral structure is present. It also informs the small-scale problem by 
ascertaining the properties (in terms of masses, magnetic fluxes, and spins) of the bound 
clouds that form, which affect the star formation efficiency. A complete theory of cloud 
formation should start from initial conditions that would be left by the previous generation 
of clouds, in which feedback from star formation has taken its toll. Clouds may be partially 
destroyed by photoevaporation and partly disaggregated by stellar winds and supernovae. 
We defer this more complex treatment for the future, here taking the first step of considering 
cloud formation in a medium that is initially relatively uniform in its large-scale properties. 

The plan of this paper is as follows: In §2, we lay out our numerical methods and the 
specifications of the models we shall study. In §3, we present results for two-dimensional 
"thick disk gravity" models, in which both magneto- Jeans instability (see Paper I) and the 
"wiggle instability" of Wada & Koda (2004) can develop, depending on parameters. In 
§4, we present the results of three-dimensional stratified models, identifying the conditions 
under which spurs and clouds form via various different self-gravitating mechanisms (but not 
through the "wiggle instability"). In §5, we conclude with a summary of our new results, 
and a discussion of their implications within the larger context of current work on galactic 
structure and star formation. 



2. Numerical Methods and Model Pcirameters 

In this paper, we study the nonlinear evolution of vertically stratified, differentially 
rotating, self-gravitating, magnetized galactic gas disks under spiral stellar potentials using 
local, three-dimensional numerical simulations. Similar work in a two-dimensional, razor- 
thin disk geometry was reported in Paper I, while three-dimensional evolution of gas in 
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the absence of the spiral potential perturbation was studied in Paper II. In this section we 
briefly summarize the numerical methods and model parameters that we adopt; the reader 
is refereed to Papers I and II for the more complete description. 



2.1. Basic Equations and Numerical Methods 

We study the nonhnear evolution of self-gravitating gas under the influence of a stellar 
spiral potential that is tightly wound with a pitch angle i <^ 1, and rigidly rotating at 

a constant pattern speed Qp with respect to an incrtial frame. For local simulations (i.e. 
domain <^ Rq) involving spiral arms, it is advantageous to construct a local Cartesian frame 
centered on a position Rq, 0o = ^pt, Zq = 0, that corotatcs with the stellar spiral pattern. 
The local Cartesian frame is inclined with respect to the R, coordinate directions by an 
angle i in such a way that x and y correspond to the directions in the midplane perpendicular 
and parallel, respectively, to the local segment of the spiral arm, while z denotes the direction 
perpendicular to the galactic plane (Roberts 1969; Balbus 1988). The simulation domain 
is a rectangular parallelepiped with size x Ly x L^. The dimensions of the box are 
L-r = 27tRo sin i/m, equal to the arm-to-arm distance for an m-armed spiral, Ly = 2Lx, and 
Lz = 8Hq, where Hq denotes the vertical scale height of the gas distribution when the spiral 
potential perturbation is absent. 

In this local frame, the background velocity arising purely from galactic rotation is 
approximately given by 

Vo = -Ro(^^o - fip) sinix + [Ro{^lo - ^p) - qo^ox]^, (1) 

where Qq is the angular velocity of gas at Rq in the inertial frame and go = — {din Q/ din R)\ro 
is the local shear parameter of the background flow in the absence of the spiral arm (Paper 
I). For a flat rotation curve, Qq = 1. We assume that the gas velocity v induced by the spiral 
potential is much smaller than RqHq, and neglect terms arising from curvature effects in the 
coordinates (e.g., Goldreich & Lynden-Bell 1965b; Juhan & Toomre 1966). In this local, 
rotating frame, the MHD equations are written as 

^ + V- (pyt) = 0, (2) 

&v 11 

— + vt • Vv = — VP + — (V X B) X B + go^^o^^xY - 212o x v - V($s + *ext), (3) 
at p 4:np 

_ = Vx(vtxB), (4) 

V'$s = 471(7^, (5) 
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and 




(6) 



(cf. Roberts 1969; Roberts & Yuan 1970; Shu, Milione, & Roberts 1973; Balbus & Cowie 
1985; Balbus 1988), where Vt = Vq + v is the total velocity in the rotating frame, Cg is 
the isothermal sound speed, $s is the self-gravitational potential of the gas, and $ext is the 
external stellar potential. Other symbols have their usual meanings. In all the simulations 
presented in this paper, for simplicity we adopt an isothermal equation of state in both space 
and time, as expressed by equation (6). 

The imposed external stellar potential $ext varies both in the plane and perpendicular 
to the plane; to lowest order it is separable into two parts as 



The first, quadratic term in z describes the variation in the gravitational potential near 
the midplane from a self-gravitating stellar distribution with surface density S^, and vertical 
velocity dispersion a^^z-^ 

Since the bulk of gas remains within one scale height of the stellar disk, this is quite a 
good approximation for the stellar potential in studying the dynamical evolution of the gas. 
The second, sinusoidal term with amplitude $sp in equation (7) represents a local analog of 
the logarithmic spiral potential used in Roberts (1969) and Shu, Milione, & Roberts (1973). 
Since $sp < and \x\ < the spiral potential attains its minimum at the center {x — 0) 

of the box. 

We integrate the time-dependent, ideal MHD equations (2)-(7) using a modified ver- 
sion of the ZEUS code (Stone & Norman 1992a,b). It uses a time-exphcit, operator-split, 
finite-difference scheme for solving the MHD equations on a staggered mesh, and employs 
the "constrained transport" and "method of characteristics" algorithms to maintain the di- 
vergence free condition of magnetic fields as well as to ensure accurate propagation of Alfven 
waves. For less diffusive transport of hydrodynamic variables, we apply a velocity decom- 
position method in updating Vy (Kim & Ostriker 2001). We also employ the Alfven limiter 
algorithm of Miller & Stone (2000), setting the limiting speed of the displacement current to 
ciim = 8cs; this value of allows a good dynamical range in the low-density, high- 2; regions. 



^The formula for the gravitational potential of a self-gravitating stellar disk treated as an isothermal 
fluid is = 2cr^^ lncosh(2/iJ*), where /f* = c7^^/(7rGS*) (cf. Gilmore, King, & van der Kruit 1990). 

The central stellar density is po,* = S*/(2iJ*). For z/H^ <C 1, ^*{z) « {at.^z/H^.Yz'^ = 2'!rGpo,*z^ = 
{nGi:,/a,,,fz^. 
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We solve the Poisson equation by combining the fast Fourier transform method in sheared 
horizontal coordinates (Gammie 2001) with the Green function method for vertical integra- 
tion (Miyama, Narita, & Hayashi 1987). We adopt the shearing box boundary conditions of 
Hawley, Gammic, & Balbus (1995) in which the x-boundaries arc shearing-periodic and the 
y-boundaries arc perfectly periodic, while implementing the outflow boundary conditions of 
Stone & Norman (1992a) at the ;2-boundaries. 

We have parallelized the numerical code we use with both Open-MP and MPl, for use on 
both shared- and distributed-memory platforms. The parallelization is achieved by domain 
decomposition along the 2;-dircction, which is convenient for our hybrid technique of solving 
the Poisson equation. Our high resolution models in three dimensions have 128 x 256 x 128 
zones in {x, y, z). 

2.2. Model Pcirameters 

The ISM in galaxies is highly inhomogcncous, turbulent, and multi-phase (e.g.. Field, 
Goldsmith, & Habing 1969; McKee & Ostrikcr 1977; Hcilcs 2001; Heiles & Troland 2005; 
Wolfire et al. 2003), and a fully realistic treatment of ISM evolution entails consideration of 
heating, cooling, and other physical processes that may significantly affect the density and 
temperature structures of the gas (e.g, Vazquez- Semadeni, Gazol, & Scalo 2000; Kritsuk 
& Norman 2004; Piontek & Ostriker 2004; Audit & Hennebelle 2005; Piontek & Ostriker 
2005). In order to focus on the role of self-gravity in forming intermediate-scale spiral-arm 
substructure, we instead consider for all models presented in this paper homogeneous initial 
gas distributions, adopting a simple isothermal equation of state with an "effective" speed 
of sound Cs- The effects of ISM heating/cooling and the associated multi-phase cloudy gas 
distribution on the interaction with spiral arms will be considered in a subsequent paper. 

Our initial disks, in the absence of the spiral perturbation (<l>sp = in eq. [7]), are 
vertically stratified with density Pq{z), have a uniform surface density Sq = / po{z)dz, and 
are threaded by magnetic fields Bq = BQ{z)y pointing parallel to the spiral arm. The 
initial equilibrium satisfies force balance between the total (thermal plus magnetic) pressure 
gradient and the total (self plus external) gravity along the z-direction. For the strength of 
the external vertical gravity, we define Sq = (a^^z'^o)^ / {cgT^*)^ and take Sq — 1- This choice 
of So corresponds to average galactic-disk conditions of a^^z ~ 20 km s~^, Ri 35 Mq pc~^, 
and total gas surface density Sq ~ (11 — 16) M© pc~^ (e.g., Kuijken & Gilmore 1989; 
Holmbcrg & Flynn 2000). Physically, Sq measures the ratio of gaseous self-gravity to stellar 
gravity at one scale height Hq of the gas; Hq = So/(2po[0]) ~ 200 pc for the adopted sets of 
parameters in this paper (see below and Paper II for a discussion of the vertical equilibrium 



-8- 



specification) . 

The three key dimensionless parameters that characterize our simulation models are 



^""TrGEn "^'^(r-Okms-O (sekms-ikpc-i) (l3 M^pc-O ' 



ttGEo V7.0kms-V V36kms-i kpc-iy \13 MqPc 

c[ ^ 47rpo(z) 
vl Bl{z) 



c2 47rpo(z)c^ 



and 



m ( |$sp| 



sini \ i?o^o 



(10) 



where f a = -B/(47rp)^/^ is the Alfven speed. The Toomre stability parameter Qo describes the 
gas surface density relative to various critical values that mark thresholds for axisymmetric 
gravitational instability. For a razor-thin disk, the instability threshold is at Qc = 1 (Toomre 
1964); for a purely self-gravitating (sq = oo), vertically-stratified disk, the threshold is at 
Qc ~ 0.676 (Goldreich & Lynden-Bell 1965a; Gammie 2001); and for a disk subject to both 
self-gravity and external gravity with sq — 1, the threshold is at Qc ~ 0.75 (Paper II; Kim, 
Ostriker, & Stone 2003, hereafter Paper III). The dimensionless parameter F in equation 
(10) measures the ratio of the maximum perturbation force from the external spiral potential 
to the mean radial gravitational force (Roberts 1969). Initially, we take (3q to be constant 
everywhere, so that Bo{z) oc [po{z)Y^'^. 

To represent solar neighborhood conditions in dimensional quantities, we use galacto- 
centric radius Rq = 10 kpc. With Qq = 26 km s~^ kpc~^, the local epicyclic frequency is 
Ko = {R-^d[R*n^]/dR)\R^y^^ = (4 - 2go)^/^f^o ~ 36 km s^^ kpc^^ for a nearly fiat rotation 
curve with go ~ 1 (Binney & Tremaine 1987). The corresponding galactic orbital period is 
torh = 27r/r2o = 2.4 x 10® yr {Qq/26 km s^^ kpc^^)"^, which we adopt as a fiducial time unit 
in our presentation. Our results can be rescaled to use other galactic parameters, provided 
that the dimensionless ratios Qp/Qo, L^/Rq, and ^qRq/cs (as well as go, Qo, Po, and F) 
remain the same. 

For spiral arm parameters, we adopt for all our models pattern speed Qp = flo/2, pitch 
angle sini = 0.1, and azimuthal wavenumber m = 2; the corresponding arm-to-arm distance 
is Lx — 2Tr Ro sin i/m. For our fiducial parameters, the quasi-radial size of the simulation 
domain is therefore — 3.14 kpc. 

To simulate various galactic conditions, we select the following three sets of the param- 
eters: {Qo, Po, F) = (1.8, oo, 5%), (1.5, 10, 5%), and (1.2, 1, 10%). For our fiducial Ro and 
Qo, and taking Cs — 7 km s~^, the corresponding gas surface density, total gas mass Mtot con- 
tained in the simulation box, and scale height Hq are 11 M© pc~^, 2.3 x 10® Mq, and 200 pc 
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for the Qq = 1.8 model with /3o = oo; 13 M© pc'^ 2.8 x 10^ M©, and 170 pc for the Qo = 1-5 
model with /5o = 10; and 16 Mq pc'^ 3.4 x 10^ M©, and 180 pc for the Qq = 1.2 model 
with f3o = 1. For the f3o = 1 and 10 models, the respective midplane magnetic field strengths 
are 4.4/xG and 1.3;uG, respectively. The model parameter sets are chosen to ensure that all 
the models are gravitationally stable to quasi- axisymmetric (wavefronts parallel to the spiral 
arm) perturbations even when the effect of the spiral potential is included. Our objective in 
this paper is to explore how nonaxisymmetric disturbances in these systems evolve, subject 
to interactions with spiral arms. We also briefly consider, in §4.2.4, development of a model 
which is unstable to quasi-axisymmetric disturbances. 



3. Two-Dimensional Models with Thick-Disk Gravity 

Paper I studied the formation of spiral-arm and interarm substructure in two-dimensional 
disks that were treated as being infinitesimally thin. One of the drawbacks of the thin-disk 
approximation of Paper I is that it overestimates self-gravity at the disk midplane for modes 
whose wavelengths approach the disk scale height (e.g., Toomre 1964; Elmegreen 1987). 
For a given strength of the spiral arm potential, therefore, razor-thin disks tend to have 
larger density enhancement at the spiral shock than vertically-resolved disks, increasing the 
susceptibility to gravitational instability. Consequently, the critical F values for stable, one- 
dimensional spiral shock configurations are lower in razor-thin disks than in more realistic 
vertically extended disks; to ensure gravitational stability to quasi-axisymmetric perturba- 
tions, razor-thin disk models considered in Paper I were limited to F < 3% for a range of 
Qo and /?o values. 

We revisit two-dimensional disk models here, now taking into account the geometrical 
dilution of self-gravity due to finite disk thickness. Instead of the standard kernel appropriate 
for a thin disk, we use a "thick-disk" gravitational kernel such that the potential and density 
Fourier modes with wavenumber k at the midplane are related by 

and ^s{k) = for k = 0, where the disk scale height Hq is held fixed in both space and 
time. Equation (11) is in fact exact for an exponential density distribution p oc e~^/^° 
(e.g., Elmegreen 1987), and generally yields a result within 15% of the exact solution of the 
three-dimensional Poisson equation for self-consistent vertical equilibria (Paper 11). For our 
two-dimensional models, the dependences of all fluid variables on the vertical coordinate are 
neglected. The results of these two-dimensional disk models allow us to quantify the impact 
of thick-disk gravity on the formation of arm/interarm substructure without the potential 
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influence of dynamical instabilities that critically involve the vertical dimension. 

We have run three sets of models with differing numerical resolutions. The model 
parameters and simulation results are given in Table 1. Column (1) labels each run. Columns 
(2) and (3) list the Toomre Qq parameter (eq. [8]) and the field strength in terms of the 
plasma parameter /3o (eq. [9]) of an initial disk when a spiral arm perturbation is absent, 
while column (4) gives the amplitude of the spiral potential in terms of F (eq. [10]). Column 
(5) indicates the numerical resolution A^^ x Ny of the model. Since our simulation box 
has Ly = 2Lx, Ny = 2Nx makes each cell square in the x-y plane. The peak surface 
density Egp of the resulting spiral shock configuration is given in column (6), while column 
(7) gives the corresponding local value Qsp = (5o(Ssp/So)^^/^ at the spiral density peak; 
this scahng follows from the constraint of potential vorticity conservation (Hunter 1964; 
Balbus 1988; Gammie 1996; Paper I). The width W of the gaseous spiral arm defined at 
E = (max[E] + min[E])/2 is listed in column (8). Finally, the mean separation (A^^) of the 
structures that develop in each model is given in column (9) in terms of the arm-to-arm 
distance and in column (10) in terms of the normalized wavcnumber Ky = Xj^sp/\, 
where the local, thin-disk Jeans wavelength at the arm density peak is defined by 



Note that structures that form from magneto- Jeans instability (MJI) in a featureless, uni- 
form, low-shear, razor-thin disk favor Ky — 0.50 — 0.75 (Kim & Ostriker 2001), while MJI- 
driven fragmentation occurring inside spiral arms in the razor-thin limit was found to have 
~ 0.45 - 0.54 (Paper I). 

To initiate our models, we begin with a disk having a uniform shear profile expressed by 
equation (1), and uniform surface density and magnetic fields corresponding to the cho- 
sen values of Qq and /3o- Using a one-dimensional grid in x, we turn on the external 
spiral potential and slowly increase its amplitude up to a desired level, F. This yields a 
one-dimensional equilibrium spiral shock profile, which we then use to initialize our two- 
dimensional simulations. On top of the background profile, we apply density perturbations 
created by a Gaussian random field with flat power for 1 < kLx/Iji < 128 and zero power 
for kLx/2TT > 128. The standard deviation of the density perturbations is flxed to be 3% 
in real space. ^ We then follow the two-dimensional flow as perturbations evolve and grow 
through interaction with the background spiral shock, eventually forming self-gravitating 



^This flat power spectrum with low amplitudes by no means represents realistic interstellar perturbations; 
it is chosen simply to allow identification and study the linear and nonlinear responses of the most dominant 
modes of the system. The ultimate outcomes of the simulations are not, however, sensitive to this choice. 
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substructures. Note that an alternative initialization procedure is simply to start with a 
uniformly-shearing disk, and slowly turn on the spiral forcing term; this is the procedure we 
follow for our three-dimensional models (see §4), and we have also tested two-dimensional 
models with this method. We find, for our two-dimensional models, essentially the same end 
results regardless of the initialization procedure. 

3.1. Magneto- Jeans Instability inside Spiral Arms 

In this subsection we describe the evolution of model ME2dl (and its lower-resolution 
counterparts ME2d2-3), which have equipartition magnetic fields. In these simulations, 
gaseous spurs and gravitationally bound clouds form as a direct consequence of MJI operating 
within spiral arms. The phrase 'bound clouds' in this paper refers to gaseous clumps that 
would collapse in a runaway fashion if physical processes such as turbulence, fragmentation, 
or star formation did not subsequently limit this collapse. We note that if the razor-thin 
gravitational kernel had been used (i.e. Hq — in eq. [11]), then the spiral shocks produced 
by this set of model parameters (Qo = 1-2, Po — 1, F — 10%) would have been gravitationally 
unstable to quasi-axisymmetric perturbations. Instead, with nonzero disk thickness, a quasi- 
axisymmetric, stable shock configuration is possible at a value of F = 10%. 

Figure 1 plots time evolution of the maximum surface density of model ME2dl (together 
with other two-dimensional models with different field strength) , while Figure 2 shows snap- 
shots of model ME2dl at t/torb = 2.5 and 3.2. The perturbations introduced into the fiow 
shear around and relax initially (t/torb ~ 0.4), and then various modes begin to grow due to 
self-gravity (t/torb ~ 2). In the absence of magnetic fields, growth of perturbations under the 
reversed-shear conditions within the arm can be inhibited by Coriolis forces; these usually 
reduce the mass-collecting effect of self-gravity in a rotating system. However, magnetic ten- 
sion forces from embedded field lines counteract Coriolis forces, so that potential vorticity 
is no longer conserved and perturbations can grow rapidly (Lynden-Bell 1966; Elmegreen 
1987; Kim & Ostriker 2001). The MJI process works best when shear is weak. The spatially 
varying sense of shear inside a spiral arm (reversed, then "normal") keeps the overall shear 
rate close to zero, therefore providing favorable conditions for the development of MJI within 
the arm (Paper I). 

Since the perturbations in model ME2dl initially have low amplitude, they remain in the 
linear regime through t/toih ~ 1-5. It is not until t/torb ~ 2, when perturbations have crossed 
the spiral arm twice (since tcross = ^orb/[?7^(l — ^p/^)]), that the most dominant MJI mode 
emerges and shapes the condensed gas fiowing into the interarm region into spur structures. 
Figure 2a shows the surface density (in logarithmic scale) and velocity vectors at t/torb = 2.5 
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viewed from a frame corotating with the spiral arm. We identify 4 spurs that protrude, at 
fairly regular intervals, perpendicularly from the main spiral arm, then sweeping back into a 
trailing configuration in the interarm region. The normalized wavenumber corresponding to 
the spur separation A^, in model ME2dl is Ky = Aj^sp/\ = 0.11. These spurs move in the y- 
direction with a speed Vy ~ O.SOi^o^o relative to the spiral arm, implying that in the inertial 
frame they follow very closely the background galaxy rotation at the arm center. Since the 
lower-resolution model ME2d2 with 128 x 256 zones also forms 4 spurs, while model ME2d3 
with 128 X 128 zones results in 3 spurs (see Table 1), we conclude that the number of spurs 
that form in our two-dimensional thick-disk models is independent of numerical resolution 
as long as the y-dimension of the simulation domain has 256 zones or more. 

Figure 2b draws selected gas streamlines (in red) seen in the stationary-spur frame as 
well as the surface density of model ME2dl at t/^orb — 2.5. In this figure, coordinates are 
transformed such that the left boundary corresponds to the initial location of the spiral 
shock front. Because the background flow is shearing and expanding/contracting, the x- 
wavenumbcr of Lagrangian perturbations that move away from the shock front varies as 
kx — —Tky with 



dimcnsionless elapsed time that is measured from the shock location (or density peak), Xgp, 
ky = 2n/\y with Aj^ corresponding to the spur spacing, and kx{0) the initial x- wavenumber 
at T = (Paper I; see also Balbus 1988). Figure 2b also plots (in black) the theoretical 
wavefronts of spurs given by dy/dx — —k^/ky — T with an initial condition Kx{Q) = 
^a;(0)/A".j sp = 0.45, where fcj sp ~ 27r/Aj,sp. The fact that the shape of spurs matches quite 
well with the theoretical prediction suggests that the former simply reflects the shearing 
and expanding properties of the background flow (Paper 1). Note that the gas streamlines 
rather quickly converge to the spur wavefront as they move downstream from the spiral 
shock, indicating that spurs grow stronger by gathering material mainly along the y-direction 
(parallel to the spiral arm). 

When the density within the spurs has grown sufficiently, self-gravity causes them to 
fragment into gravitationally bound condensations. Figure 2c shows the density and mag- 
netic field lines at t/torb = 3.2 of model ME2dl. The magnetic fields roughly parallel the 
arm overall, although they pinch inward within the spurs and are strongly twisted locally in 
the vicinity of the bound clumps. Bending of field lines is most severe in interarm regions 
where the gas moves faster. Model ME2dl forms 4 clumps with mass M ~ 3.3 x 10''' Mq 
each; this is about an order of magnitude larger than the clumps that formed in razor-thin 
disk models of Paper I. Roughly 40% of the total mass, therefore, is collected into bound 
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where TZ 



Sgp/E is the local surface density expansion factor, r = v^\:dx is a 
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clumps. These condensations are magnetically supercritical, with the mean mass-to-flux ra- 
tio M/$B ~ 2.0G-^/^ where $b is the magnetic flux that passes through each clump and G 
is the gravitational constant.^ The numerical box of model ME2dl initially has a mass-to- 
flux ratio of M/^b = S.6G~^^^, about twice as supercritical as the clumps that form. This 
result is consistent with Vazquez-Semadeni et al. (2005a), who found that self-gravitating 
substructures formed in turbulent MHD simulations are in general less supercritical than the 
parent system, presumably because fragmentation occurring along the flux tubes reduces the 
cloud mass while preserving magnetic flux (see also Li et al. 2004). Various physical prop- 
erties of bound clouds that form in two-dimensional thick-disk models are quite similar to 
those in full, three-dimensional models; we defer detailed discussion to §4.2.3, where the 
vertical stratiflcation of disks is explicitly taken into account. 

The reduced self-gravity in disks with flnite thickness results in a smaller number of 
spurs compared to a razor-thin disk with properties otherwise the same. In order to check 
whether our simulation results are consistent with the linear theory prediction, we perform a 
linear stability analysis in a Lagrangian frame comoving with the background flow through a 
spiral arm. Assuming that the perturbed quantities are well described by plane waves with 
sinusoidal variations on scales ^ R, one can show that equations (2)-(6) and (11) lead to 
the following set of linearized equations: 
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where Sa = 5S/S, Su = iSvxkj^sp/^o, Sv = iSvykj^sp/^o, Sm is the normalized perturbed 
vector potential, Ky = ky/kj^sp, K = |i^^y|(l -l- T'^Y^'^ is the total wavenumber, and a = 
{cskj^spl^oY (see Paper 1). 

We choose Hq = 180 pc for thick-disk gravity and adopt the equilibrium density and 
velocity profiles (in x) of model ME2dl as a background state. By taking Sa = 1 and 
Su = Sv = Sm = as an initial condition and varying Ky and kx{0)/ky, we integrate 
the perturbed equations in time. The resulting amplification magnitude T = log |5Emax/S| 



■^The critical mass-to-flux limit is 0.16G 
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and the growth time tgrow at which maximum amphfication occurs are plotted in Figure 3 
as sohd and dotted contours, respectively. Figure 3 also marks as a rectangular box the 
parameters (-ft'i^p] = 0.45 and Ky = 0.11) that give the best fit to the shape and spacing 
of spurs formed in the simulation of model ME2dl (the size of the rectangle represents the 
uncertainty determining Kx{0) or Ky). 

Figure 3 shows that Ky = 0.11 is in good agreement with the expectations from linear 
theory for the wavelength parallel to the spiral arm of the dominant mode. However, the 
spurs that dominate the simulation in model ME2dl have an x-wavenumber at the spiral 
shock front that is a factor 2.5 times larger than that which would yield the largest am- 
plification, based on a linear-theory analysis. This is because the initial perturbations in 
model ME2dl have such low amplitudes that two traversals of an arm are required before 
the perturbations grow into the nonlinear regime. Over this time, K^ of any initial perturba- 
tion varies kinematically due to compression, expansion, and shear of the background flow, 
and overall K^ increases due to the large shear in the interarm region. Perturbations that 
are amplified (via MJl) during a first passage through a spiral arm become preferentially 
trailing with high Kx{0) before they enter the next spiral shock (see Figure 2a). Because 
of the strong compression immediately behind the shock, -^'^^(O) of these trailing wavelets is 
further increased before they enter the next stage of amplification by MJL Consequently, the 
perturbations that eventually grow to emerge as spurs in model ME2dl have quite a large 
value oi Kx:{0). 

To quantify how finite disk thickness is expected to affect the development of MJl 
inside spiral arms, we vary Hq and search for i^a;,max(0) and Ky^^^.^ that give - based on 
linear theory - the maximum amplification magnitude Fmax for a given value of Hq. Figure 
4 plots the resulting i^^j;,max(0) and Ky^^^ as solid lines as well as Fmax as a dotted line. 
The equilibrium surface-density profile of model ME2dl (but with varying Hq) is again 
used. The hnear theory suggests that the MJI inside spiral arms in a disk with Hq — 180 
pc requires i^^y,max to be 2.6 times smaller than in a disk with Hq — 0. The corresponding 
amplification magnitude Fmax is 2.3 times smaller for thick disks compared to thin disks. For 
comparison, we note that in the razor-thin models of Paper 1, spurs identified in simulations 
with /3o = 1 typically have Ky ~ 0.45 — 0.54 (or ~ 1.8 — 2.2Aj^sp); with larger Ky (or 
smaller A^^) corresponding to smaller-Qo models. Since the amplification factor is enormous 
when Ho — 0, the modes that come to dominate in razor-thin disks tend to be selected by 
compromise between maximum amplification and earliest growth (Paper I), and thus have 
Ky larger than that corresponding to i^^y,max- 

In the thick-disk models with their lower amplification, on the other hand, early growth 
gives httle advantage over other modes. Model ME2dl therefore chooses the mode corre- 
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spending to the largest total amplification, whose wavelength turns out to be Xy ~ 9Aj,sp 
which is ~ 4 — 5 times larger (in terms of Aj^/Aj_sp) than that found for razor-thin disk 
models in Paper I. Models with thick-disk gravity have peak surface densities twice as large 
as their razor-thin counterparts with the same F and Qq, which reduces the Jeans length 
by about a factor of two. Thus, for cloud mass scaling as A'^Sgp oc (Aj^/Aj,sp)^S^^, one can 
expect that bound clouds resulting from spur fragmentation in thick-disk gravity models 
should be ~ 8 — 10 times more massive than in razor-thin disks with the same F and Qq, or 
~ (3 — 4) X 10^ Mq in dimensional units. This is entirely consistent with the results of our 
numerical simulations. 



3.2. The Wiggle Instability 

We now turn to model MS2dl (and its lower-resolution counterparts MS2d2-3), with 
Qo = 1.5, F = 5%, and (3o = 10 (i.e., va = Cs/vTo), and unmagnetized model H2dl (with 
counterparts H2d2-3) having Qo — 1.8, F — 5%, and Po — oo. As Figure 1 shows, in models 
MS2dl and H2dl perturbations immediately begin to grow almost exponentially over time, 
while strong growth in model ME2dl does not begin until after two successive passages 
through spiral arms. The growth rates of the maximum surface density in models MS2dl 
and H2dl arc almost identical, ~ 0.70i7o) suggesting that moderate magnetic fields are not 
important at least in the linear stages of growth in these simulations. 

The two models MS2dl and H2dl have similar evolution, but this is quite distinct 
from the behavior of equipartition- magnetization model ME2dl. Unlike model ME2dl, in 
which the spiral shock remains relatively straight during the initial phase of MJl growth, 
models MS2dl and H2dl show small-scale instability that wiggles the spiral shock front and 
forms structures with high vorticity. The left panel of Figure 5 plots the surface density 
map in logarithmic gray-scale of model MS2dl at t/torh — 1-6. The dotted line marks 
the location of the spiral shock front, while the dashed line represents the position where 
shear in the background flow vanishes {dvT,y/dx — 0). Evidently, the vortical clumps that 
emerge are closely connected to the shock discontinuity. Column (9) in Table 1 shows that 
the hydrodynamic models H2dl-3 form more clumps than the weakly magnetized models 
MS2dl-3, and that the number of vortical clumps produced depends rather sensitively on 
numerical resolution. 

At t/torh ~ 1.7, the growth of vortical structures in model MS2dl saturates due to mag- 
netic tension from bent field lines as well as nonlinear effects. (The flow properties inherent 
in the background spiral shock may also prevent the clumps from growing further.) At this 
time, the clumps are not sufficiently self-gravitating to experience gravitational collapse due 
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to their own weight. The clumps wander shghtly inside the postshock regions and coUide 
with each other to form bigger clumps. As Figure 56 shows, mergers of 10 low-density clumps 
at t/torb = 1-6 result in 5 high-density, gravitationally bound clumps at t/torb = 2.7. These 
bound clumps have an average mass ~ 1 x 10^ M©, a factor of 3 smaller than that origi- 
nating from MJI in model ME2dl, and are supercritical with an average mass-to- flux ratio 
M/$B ~ 3^-^/2 xhe interarm spur structures prominent in model ME2dl (see Fig. 2) are 
not clearly visible in models MS2dl and H2dl. 

What causes the spiral shocks to wiggle and form vortical structures in models MS2dl 
and H2dl? Figure 6 shows the temporal evolution of gas surface density and magnetic 
field lines in the rectangular section shown in the left frame of Figure 5. Note that the 
vortex generation and evolution of the magnetic field topology in model MS2dl are similar 
to those that result from Kelvin- Helmholtz instabihties in unmagnetized (e.g., Corcos & 
Sherman 1984) or magnetized (e.g., Malagoh, Bodo, & Rosner 1996; Prank et al. 1996) 
shear layers. This suggests that the vortical clumps in model MS2dl may arise from Kelvin- 
Helmholtz instability in the sheared flow - i.e. azimuthal streaming - associated with spiral 
shocks. Wada & Koda (2004) performed non-self-gravitating global simulations of razor-thin, 
unmagnetized disks and found that when spiral perturbations are strong and have large pitch 
angles, the shock front that develops wiggles and forms discrete clumps bearing remarkable 
resemblance to those found in our models. Wada and Koda termed this process a "wiggle 
instability". Based on the low Richardson numbers in the postshock regions, they argued 
that the formation of these clumps could be due to Kelvin-Helmholtz instability. As Wada 
and Koda noted, however, the Richardson number criterion is only a necessary condition for 
stability (Chandrasekhar 1961), and should not be regarded as an instability criterion. The 
critical value of Ri may, in addition, be affected by rotation. Thus, while it is likely that 
the wiggle instability is indeed a manifestation of Kelvin-Helmholtz instability in the shear 
flows produced by spiral shocks, it is not yet certain. 

Since models considered in Wada & Koda (2004) are non- self-gravitating and unmag- 
netized, the wiggle instability, whether it is related to Kelvin-Helmholtz instability or not, 
relies neither upon magnetic fields nor self-gravity. Paper 1 found that gas in razor-thin disks 
with F < 3% is subject to either the MJI or the swing amplifier, but is stable to the wiggle 
instability; i.e., non-self-gravitating models at these values of F are expected to be stable. 
We have run a number of two-dimensional simulations of non-self-gravitating, unmagnetized 
disks with varying F (not listed in Table 1), and found that the wiggle instability becomes 
manifest when F > 4%. This suggests that spiral shocks must be fairly strong in order 
to trigger the wiggle instability^. Since non-self-gravitating gas behind a spiral shock is 



Wada & Koda (2004) in fact set up extremely strong spiral shocks whose amplitudes amount to F « 
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locally Kelvin-Helmholtz stable due to the expanding radial velocity inside the spiral arm 
(Dwarkadas & Balbus 1996), the wiggle instability is most likely a consequence of interaction 
between preshock gas and compressed postshock flows that stream at large relative veloc- 
ity parallel to the shock front. It is also possible that spiral shocks are so strong that the 
postshock velocity perpendicular to the arm is vanishingly small, limiting the stabilization 
to Kelvin-Helmholtz modes. In any event, we do not believe the wiggle instability is likely 
to be important in spiral arms of real galaxies, because as we will show in the next section 
it does not occur in three-dimensional models, in which the fluid variables are not constant 
with height. 



4. Three-dimensional Models 

To investigate the effects of disk stratiflcation and other vertical variations on the forma- 
tion of spiral-arm substructures, we have performed three sets of three-dimensional numerical 
simulations. The chosen parameters for the three-dimensional models listed in Table 2 are 
identical to those of the two-dimensional models presented in the previous section. We flrst 
set up an axisymmetric disk with equilibrium density distribution po{z) consistent with a 
chosen set of Qo, Po, ^-nd Sq. We then apply initial perturbations to the background den- 
sity po{z), using a spatially uncorrelated, Gaussian random field that has fiat power for 
1 < kLx/2TT < 64 and zero power for 64 < kLx/27c. For the amplitude of the perturbations, 
we fix the standard deviation of perturbed density in real space to be 3% of the initial mid- 
plane density. Next, we slowly turn on a spiral potential perturbation (the second term in 
eq. [7]) such that it acquires the full strength F at t/torh = 1-5. 



4.1. Vertical Structure of Spiral Shocks 

Figure 7 shows evolution of the maximum density in high-resolution three-dimensional 
models ME3dl, MSSdl, and H3dl. As the strength of the spiral potential increases, a 
density wave grows rapidly in the gas, until the initially-sinusoidal profile steepens into 
a shock. During this phase, the growth of small-scale nonaxisymmetric perturbations is 
weak compared to the large-scale nonlinear response of the gas disk to the imposed spiral 
potential perturbation. At t/tovh ~ 1.6 — 1.8, slightly after the spiral potential attains its full 
strength, the spiral shock reaches its maximum strength, leaving density and velocity fields 



2ecoti ~ 110% for the parameters e = 0.1 and i = 10° in their model A; the arm/interarm contrasts reached 
two orders of magnitude. 
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that fluctuate around the respective average conflguration. At this stage, the spiral shock 
still remains almost axisymmetric. 

Paper I and §3 of the present paper showed that when vertical degrees of freedom 
are ignored, quasi-axisymmetric, steady-state profiles of gaseous spiral shocks can be read- 
ily obtained through time-dependent calculations in an external spiral potential (see also 
Woodward 1975). This implies that onc-dimcnsional spiral shock profiles in the galactic 
midplane represent stable shock equilibria (provided that gaseous self-gravity and external 
spiral forcing are not too strong; see Paper I). When we allow for vertical degrees of free- 
dom in the fluid variables, however, we find that the profiles that develop are not in general 
quasi-steady. In the initial stages of all the three-dimensional models studied in this work, 
quasi-axisymmetric, two-dimensional spiral shocks develop in the x-z plane (hereafter XZ 
spiral shocks) that are in general non-stationary, swaying loosely back and forth in the 
direction perpendicular to the arm. 

Based on our simulation results, the "fiapping" motions of XZ spiral shocks are strongest 
at high altitudes, and stronger in model ME3dl than in models MSSdl and H3dl. We 
measure the maximum shock front excursion from its mean location as about O.O6L3;, 0.02Lx, 
and O.OSLj. a,t z/Hq — 1 in models MESdl, MS3dl, and HSdl, respectively. The non-steady 
nature of XZ spiral shocks that we find is in fact commonly seen in other numerical models 
(e.g., Martos & Cox 1998; Gomez & Cox 2004). Although Soukup & Yuan (1981) constructed 
stationary shock profiles in a vertically extended disk, their models neglected the effects of 
the Coriolis force, self-gravity, and magnetic fields. The fiapping motions of XZ spiral shocks 
may just be overshooting of gas caused by instantaneous force imbalance across the shock, 
or may represent a dynamical instability. Detailed discussion on the physical origin of the 
spiral shock flapping, and the characteristics of the associated velocity and density flelds, 
will be presented in a subsequent paper. 

Since spiral shocks in three dimensions exhibit temporal fluctuations, it is useful to take 
space-and-time averages in order to visualize the shock structures in the x-z plane. Figure 8 
illustrates this using model MS3dl, with the density and velocity flelds averaged along the 
y-direction and also over an orbital period from t/torh — 1-7 to 2.7. Figure 8a shows gas 
density in logarithmic scale and selected streamlines. The streamlines run almost parallel 
to the galactic midplane in the interarm region, but sharply plunge toward the midplane in 
the arm. The location of the averaged shock front is also indicated by a heavy line. The 
low vertical velocities of the gas in the interarm regions implies that vertical force balance 
is fairly well maintained there. As the material enters the spiral arm region, it is shocked 
and compressed. The increase of gas density due to the shock compression is largest at the 
midplane, which in turn produces strong vertical gravity and pulls high-altitude gas toward 
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the midplane. The vertically-plunging gas further increases the midplane density after the 
shock, setting up repulsive pressure gradients that cause it to rebound back to high-altitude 
regions. The furrow in the streamlines behind the shock front reflects this vertical dive and 
bounce. 

Figure Sb,c plots as solid lines the surface density profile E(x) = J^^p{x, z)dz and 
the corresponding gas scale height H{x) = Il{x) / {2 p[x , 0]) obtained from the averaged XZ 
spiral shock shown in Figure 8a. Also plotted as dotted lines are the and H{x) values 
of the thick-disk two-dimensional model MS2d2 that has exactly the same parameters and 
numerical resolution as model MS3dl, but fixed Hq = 170 pc. Although the surface density 
in model MS2d2 rises slightly more steeply at the shock than in model MS3dl, they are 
overall in excellent agreement. This proves that the gravitational kernel given by equation 
(11) accurately handles the dilution of self-gravity in an extended disk. 

As Figure 8 c shows, the gas within spiral arms is more compressed toward the midplane 
than in interarm regions. The tendency for gas to compress vertically in the arms increases 
as the magnetic field strength decreases. Column (8) of Table 2 indicates that the ratio 
of arm to interarm scale heights is about 0.5 for models with Pq — oo or 10, and 0.8 for 
Po — 1 models. This is because models with weaker magnetic fields have higher arm density 
enhancement, and thus relatively stronger self-gravity that vertically compresses spiral arm 
gas. We discuss the relationship of this finding to observations in §5.2. 

Finally, we remark that our models do not exhibit the "hydraulic jump" behavior that 
Martos & Cox (1998) suggested may develop, under certain circumstances, at spiral shock 
fronts. While hydraulic jumps (which yield an increase in the gas scale height where it 
is compressed in arms) can occur when the equation of state is fairly stiff, Martos & Cox 
(1998) note that when 7=1 there are only shocks. Magnetic fields parallel to the arm 
may provide extra stiffness, but self-gravity tends to draw gas toward the midplane in the 
post-shock region; in our (isothermal) models this results in a decrease of H within arms. 
It is not presently known whether the tendency of self-gravity to reduce H would overcome 
the tendency of a sufficiently stiff equation of state to increase H in arms, since the studies 
of Martos & Cox (1998) did not include self-gravity. 

4.2. Nonaxisymmetric Evolution of Three-Dimensional Disks 

4-2.1. Absence of the Wiggle Instability in Three-Dimensional Disks 

Based on our two-dimensional thick-disk models, we found that perturbations in spiral 
shocks can grow, either as a result of the MJI or the wiggle instability. One may naively 
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expect that three-dimensional models would also be prone either to MJI or to wiggle in- 
stability, not to mention to other potential three-dimensional instabilities. Surprisingly, 
however, we find that the wiggle instability seen in some of our two-dimensional models is 
absent in three-dimensional models, while MJI still operates in three-dimensional models. 
As we will show below, our three-dimensional model MS3dl with = 10 is subject to MJI, 
whereas the corresponding two-dimensional model MS2d2 (with identical in-plane resolu- 
tion to model MS3dl) was found to be unstable to the wiggle instabihty. In addition, the 
unmagnetized two-dimensional models H2dl-H2d3 developed vortical clumps via the wiggle 
instability that later merged and collapsed, but as Figure 7 shows, the corresponding three- 
dimensional model HSdl does not develop strongly overdense clumps over the entire run. 
The peak density in model HSdl fluctuates with period ~ 0.29torb, comparable to the period 
associated with epicyclic motion, ~ (2Esp/Eo)~^/^iorb ~ O.Sl^orb) within the arm peak. 

The wiggle instability seen in two-dimensional disks thus appears to be stabilized in 
three-dimensional disks, possibly due to a combination of strong vertical shear and non- 
steady flow. Figure 9 shows the vertical distribution and gradients of the azimuthal- and 
time-averaged spiral shock velocities at x/L^ = —0.05 (upstream), 0.05 (at the midplane 
density peak), and 0.22 (far downstream), from model MS3dl. The velocity shear along the 
vertical direction amounts to Qxz ^ 1 and qy^ ~ 5 at \z/H\ ~ 0.5 and Qxz ~ 2 at \z/H\ ~ 1—2, 
where qxz = ^o^\dvx/dz\ and where Qyz = ^lQ'^\dvy/dz\ . These shear rates can be compared 
with qo — 1, the radial shear rate in the background azimuthal flow, and with the maximum 
radial shear in the quasi- azimuthal flow within the arm, q — |2 — (2 — go)Ssp/So| ~ 4.3. 
Perturbations in the spiral shock of model MSSdl might try to develop nonaxisymmetric 
vortical flows similarly to those produced by the wiggle instability in two-dimensional thick- 
disk models. However, these vortical structures could not remain coherent in z because of 
the strong vertical shear in horizontal velocities; vertical motions would then mix dissimilar 
vorticity, so that the vortices would be lose their integrity. In addition, non-steady flapping 
of the shock front would also help prevent coherent structures from forming. Overall, the 
various three-dimensional effects combine to suppress nonlinear development of the wiggle 
instability. 

Since our three-dimensional models explore cases with F = 5 — 10%, our results are 
not directly comparable to cases in which spiral shocks are as strong as in Wada & Koda 
(2004) 's models. Strong spiral shocks in two dimensions are more favorable for promoting 
wiggle instability, but in three dimensions they would also result in stronger shock flapping 
motions and stronger vertical shear of the in-plane velocities, which tend to suppress the 
wiggle instability. Thus, it would require direct three-dimensional simulations with very 
large F to determine whether increasing shock strength overall stabilizes or destabilizes 
spiral shocks to "wiggle" modes. However, we consider this question of limited practical 
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importance, since real astronomical systems would in general become self-gravitating if very 
strong shocks were present. In this sense, the non-self-gravitating treatment of Wada & 
Koda (2004), with F ~ 110%, appears unrealistic given that the spiral shocks in their 
models have density jumps of a factor 100. Contrasts in real galaxies are far smaller; e.g. 
the notoriously strong spiral structure in M51 has arm/interarm contrasts of 5-6 (Garcia- 
Burillo et al. 1993; see also Rix & Zaritsky 1995; Patsis, Heraudeau, & Grosb0l 2001). This 
is because equilibrium shock profiles are only possible at moderate values of F (Paper 1). For 
instance, when Qq — 1.8 and (3 — oo, quasi- axisymmetric shock equilibria with thick-disk 
gravity exist only if F ^ 20%, beyond which shock profiles are highly transient or suffer from 
immediate quasi-axisymmetric gravitational instability (see section 4.2.4). 



4-2.2. Swing Amplification 

When shear is strong (with q ^ 0.3) and magnetic fields are absent, nonaxisymmetric 
modes are not subject to MJI but can grow via swing amplification (Goldreich & Lynden- 
Bell 1965b; Toomre 1981; Kim & Ostriker 2001). The growth of wave amplitudes via swing 
amplification is moderate unless the Qo value of the background medium is sufficiently small. 
Kim & Ostriker (2001) showed that swing amplification in a razor-thin, featureless disk 
with Qo ~ 1-3 yields low density fluctuations and cannot form gravitationally bound clumps. 
Paper 1 further found that the swing amplifler produces substructure growth in spiral arms 
only if the background spiral perturbation is relatively weak (with F < 1% in a razor-thin 
disk). In this case, the density enhancement in the shock is quite low and all the gas has 
normal shear, as opposed to the locally- reversed shear which occurs wherever S/Eq > 2. 

With F = 5%, model II3dl in the present work has strong shear reversal corresponding 
to Qmin = —4 inside a spiral arm. In addition, a relatively large value of Qo = 1.8 as well 
vertically-diluted self-gravity make swing amplification in model H3dl quite ineffective. The 
resulting nonaxisymmetric structures in the surface density of model H3dl at the end of 
the run {t/torh — 6.4) have the maximum power in the Xy — Ly/3 mode, which is only 
~ 0.26% relative to the axisymmetric mode; the associated maximum surface density is 1.46 
times the mean surface density. We thus conclude that the swing amplification mechanism 
is unlikely to prompt the formation of substructure within spiral arms in real disk galaxies 
where spiral arms are not so weak. Of course, as the models of Paper I showed, even 
without magnetic fields, if Qo is small enough, and/or if F is large enough, spiral arms can 
become unstable to fragmentation. This fragmentation in general first develops via quasi- 
axisymmetric instability, so the process is physically distinct from direct swing amphfication. 
We will discuss this alternative fragmentation process in §4.2.4. 
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4-2.3. Magneto- Jeans Instability in Three- Dimensional Disks 

In the preceding subsections, we have shown that the wiggle instabihty disappears in 
three-dimensional extended disk models, and that swing amplification is only moderate in 
unmagnetized, three-dimensional disks that have sizable spiral arm density contrasts but 
nevertheless are stable to quasi- axisymmetric self-gravitating perturbations. When a system 
is magnetized and has weak net shear, however, it is expected that perturbations should be 
able to grow via MJI. 

Figure 7 shows that the density in our magnetized models ME3dl and MS3dl indeed 
undergoes exponential growth, with the simulations eventually producing gravitationally 
bound condensations. As we will discuss below, typical wavelengths of the most unstable 
modes in models ME3dl and MSSdl are very similar to that of the MJI in two-dimensional, 
thick disks; these wavelengths are much larger than those of the wiggle instability in two- 
dimensional models (see Tables 1 and 2)^. In addition, spur morphologies and the physical 
properties of bound clouds that form from fragmentation of the spurs in three-dimensional 
models are similar to those due to MJI in two-dimensional thick-disk models. This supports 
the case that structure forms in our three-dimensional models as a result of MJI. Since 
gravity is a long-range force and insensitive to density structure at small scales, MJI (unlike 
the wiggle instability) can grow in spite of the unsteady flows of the background XZ spiral 
shock. The vertical velocity shear and flapping may delay the onset of bound cloud formation 
to some extent, but it is MJI that drives three-dimensional systems with magnetic fields into 
eventual runaway. 

As in two-dimensional models, the background fiow kinematics in our three-dimensional 
models sculpt the growing perturbations into spurs that branch out nearly perpendicularly 
from the main spiral arm and become trailing in interarm regions. Figure 10 displays surface 
density snapshots at t/torb = 5.6, 6.0, and 6.3 of model MS3dl with Pq = 10, while Figure 11 
is for model ME3dl with /5o = 1 at t/torb = 3.0, 3.1, and 3.2. The projected spur morpholo- 
gies in Figure 10a are quite similar to those of model ME2dl seen in Figure 2a, although 
the non-stationary nature of XZ spiral shocks tends to wash out the structure in three- 
dimensional model ME3dl compared to the two-dimensional model. The three-dimensional 
magnetized models all form four spurs, separated by Xy — Lx/2 — 1.6 kpc, exactly the same 



^As Table 2 indicates, the separation of spurs in our three-dimensional models depends on the numerical 
resolution. Models ME3d2 and MS3d2 with Ny = 128 form 3 spurs each, while the respective higher 
resolution model with A^,^ = 256 forms 4 spurs. Based on the results of two-dimensional models (see Table 
1), which show the same number of spurs when Ny — 256 or 512, we expect (but have not proved) that our 
three-dimensional models are converged. 
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as in model ME2dl. Thus, three-dimensional effects do not significantly alter the develop- 
ment of MJI. This also implies that two-dimensional models with thick-disk gravity represent 
an accurate and efficient means to explore the MJI without the computational expense of 
an independent vertical degree of freedom. 

To help visualize the spur structures in three dimensions, we show in Figure 12a a 

perspective volume rendering of an isodensity surface and a few selected magnetic field lines 
associated with the region marked by the rectangle in Figure 10a (from model MSSdl at 
^/^orb = 5.6 ). Shown also in Figure 12a are the density slices in logarithmic color scale and 
the velocity vectors at the midplane {z = 0) as well as at the y = 0A5Ly plane; those at the 
left edge {x = 0.5Lx) of the box are displayed in Figure 12b. Clearly, the dense part of the 
spur sticks out perpendicularly from the main spiral arm. The spur has an aspect ratio of 
1 : 0.35 : 0.16 in (x, y, z) and a total mass of ~ 7 x 10^ Mq; a similar amount of gas is in the 
part of the spiral arm shown in Figure 12a. Magnetic fields are overall parallel to the spiral 
arm, although the variation in the expanding velocity field off the arm causes some degree 
of radial excursions of magnetic fields. This results in magnetic pressure twice as large in 
the spur compared to inter-spur regions. As Figure 12b shows, magnetic field lines in the 
azimuthal-vertical plane are not strongly bent inward at high z, indicating that material is 
accumulated mainly along the (azimuthal) ^/-direction. 

As spurs develop further, they concentrate sufficiently to trigger gravitational fragmen- 
tation, yielding a few bound condensations. As Figures 10 and 11 show, fragmentation 
occurs within spiral arms as well as in interarm regions; fragments forming inside the spiral 
arms linger there, whereas those forming off the spiral arms are carried into the interarm 
regions and follow galactic rotation. At the end of the runs, model MS3dl forms 4 bound 
clouds with an average mass M ~ 1.2 x 10^ Mq, which is respectively about 6 and 1.2 times 
larger than the local, thin-disk and thick-disk Jeans masses at the arm density peak.^ These 
masses are defined by 



and 




^When the thick-disk gravity (eq. [11]) is used, one can show that the Jeans length in a disk with surface 
density S and scale height H becomes A*'^''^'^ = (27rc2iJ/GS)i/2 = {HX^y/^ when 27ri?/Aj > 1. The thick- 
disk Jeans mass applied to the peak arm density is then Mf^^^^ = Ssp(Aj'^"^^)^ = 27rc^i?sp/G, where Hgp is 
the scale height of gas at the spiral arm peak. 
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respectively. On the other hand, model ME3dl produces 3 bound clouds, one fewer than in 
the corresponding two-dimensional model ME2dl with equal magnetic field strength. This 
is because model MESdl with /9o = 1 suffers from strong XZ flapping of the spiral shock that 
tends to work against gravitational mass accumulation. Bound clouds that form in model 
MESdl have an average mass M ~ 3.0 x 10'' M©, comparable to the result from model 
ME2dl, which is about 12 times the corresponding thin-disk Jeans mass at the density 
peak. Roughly speaking, therefore, MJI in three-dimensional, magnetized spiral shocks 
turns ~ 20 — 30% of the total gas into dense bound clouds in a given epoch of gravitational 
condensation. 

Note that typical masses of bound clouds that form in the current three-dimensional 
models are very similar to those in two-dimensional thick-disk spiral arm models discussed 
in §3.1, as well as those in other three-dimensional models without spiral arms (Paper II; 
Paper III). In addition, bound clouds from different models share similar magnetic and 
rotational properties. For example, self-gravitating clouds in our three-dimensional models 
are all magnetically supercritical (by a factor 10 or more) with an average mass-to-flux 
ratio of M/$b ~ 2.4^-^/2 ^lodeX MS3dl and ~ LIG'^/^ ^^Q^gj MESdl. (Mass- 
to-flux of < O.IGG^^/^ would be subcritical.) As in two-dimensional cases, these bound 
clouds in three dimensions are less supercritical than the initial simulation boxes, which 
have M/^B = 7.0^-^2 for model MSSdl and 2.4^-^/^ foj. ^^^^^ MESdl, consistent with 
the result of Vazquez-Semadeni et al. (2005a). These results are remarkably close to that 
found in two-dimensional model ME2dl in §3.1 and those resulting from MJI (Paper II) or 
via swing amplification (Paper III) in three-dimensional disks without spiral features. 

Because of efficient magnetic torques exerted by field lines linking a cloud with the 
surrounding medium (e.g., GiUis, Mestel, & Paris 1974; Mouschovias & Paleologou 1979), 
bound clouds in models MSSdl and MESdl lose a substantial amount of angular momentum 
and thus rotate relatively slowly. The mean specific angular momentum of bound clouds 
at the end of the runs is Jz ~ 0.2Jgai for both models MSSdl and MESdl. Here, J = 
j p{r X \)d^r / J pSr with the position r and velocity v are measured relative to the cloud 
center, and Jgai = ^[^{M /HqY /12 is the specific angular momentum contained in a square 
patch of the galactic ISM. A cloud formed from that patch would preserve the value = Jgai 
if no angular momentum were lost during its formation (see Paper III). Similar specific 
angular momenta of clouds have been obtained in our other three-dimensional magnetized 
simulations (Paper II, Paper III), where we also showed that condensations in unmagnetized 
models do not lose angular momentum as they evolve, while those in magnetized models do. 

Finally, we remark on the effects of Parker (1966) instability on the formation of spiral- 
arm substructure in our three-dimensional models. The spiral potential compresses the 
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density and magnetic fields, decreasing the average f3 value inside the spiral shock from 10 
to 2.5 for model MS3dl and from 1 to 0.4 for model ME3dl. Taking /3 ~ 1 and the radial 
wavelength A^^ equal to the arm width, and neglecting the effect of cosmic ray pressure, the 
local dispersion relation of the Parker instability in a rotating disk presented by Shu (1974) 
yields ky ~ 0.5/ H for the most unstable Parker mode (this is insensitive to the specific choices 
of (3 and A^;). The corresponding wavelength is Xy = 27i/ky = AtiH pa 1.5 kpc ~ Ly/i. This 
agrees with the separations identified for spurs in models MS3dl and MESdl. However, 
because this wavelength also coincides with the scale predicted to have the largest growth 
from MJI (see §3.1 and Fig. 3), it is not immediately obvious whether Parker instability is 
important in prompting this growth. 

To investigate this further, we have examined the detailed distributions of magnetic 
fields and velocities during the period when overdense condensations begin to emerge within 
the arms. Figure 13 shows in logarithmic color scale the density in a YZ slice a,t x — 0.05Lx, 
overlaid with a vector field of perturbed velocity (relative to the mean value in the slice), and 
magnetic field lines. Evidently, there is no sign of the characteristic correlation of magnetic 
valleys/hills with overdense/underdense regions that Parker modes would produce, nor is 
there any sign of sinusoidal variations in the vertical velocities with azimuth. In part, the 
ability of Parker modes to grow may be suppressed by the strong vertical gradients in and 
Vy (as seen in Fig. 13 and also Fig. 9). Given the curved nature of the spiral shock front in 
the XZ plane, high-altitude regions can have large horizontal velocities with respect to the 
midplane gas below them, such that magnetic fields that begin to buckle vertically cannot 
maintain the shape required for Parker modes to develop. 

Thus, we conclude that magnetic buoyancy effects are probably not of major importance 
in GMC formation within spiral arms. Even without the suppression of Parker modes by 
vertical shear that we see here, previous work has shown that Parker instability alone is 
unable to form high-density clouds in regions away from arms (e.g., Kim et al. 2000; Santillan 
et al. 2000; Paper II), a consequence of it being a self-limiting process stabilized by tension 
forces from bent field lines (e.g., Mouschovias 1974). While Parker instability may help 
seed structure in early stages of gravitational instabilities, and may be essential in removing 
excess magnetic fiux from the disk, it does not appear crucial for massive cloud formation 
and subsequent star formation. 

4-2.4- Fragmentation in Strongly Unstable Arms 

We briefly discuss evolution of unmagnetized model HU3dl, with Qo = 1.5, F = 5%, 
and = oo. In this model, Qq is small enough to induce quasi-axisymmetric gravitational 
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instability of a spiral shock. As Figure 7 shows, the maximum density in model HU3dl at 
t/torh = 1-5 is larger by a factor of 3 than that in stable model H3dl with Qq = 1.8. Soon after 
the spiral perturbation attains its full strength, the system is dominated by the exponential 
growth of quasi-axisymmetric modes in the spiral shock. At around t/torb = 2.2, the collaps- 
ing post-shock layer reaches sufficient density that it begins to undergo non-axisymmetric 
gravitational fragmentation. This results in the formation of 13 bound condensations aligned 
along the shock front, with mass ~ 1 x 10^ Mq each. Note that while magnetized models 
MS3dl and ME3dl first develop spur structures that subsequently fragment into bound 
clouds in both arm and interarm regions, the bound clouds in unmagnetized model HU3dl 
do not require prior formation of spurs. These clouds form directly in the highest-density 
ridge of the spiral shock. The spacing of bound clouds in model HU3dl is A/Aj,o ^ 0.5, where 
Aj,o = cI/{GTjq) is the Jeans length in the initial featureless disk. This can be compared 
with A/Aj^o ~ 1-7 for the spur separation in magnetized model MS3dl with the same Qq 
and F parameters as model HU3dl (see Table 2)7 This demonstrates that the presence of 
magnetic fields makes significant differences in the dynamical evolution of gas in spiral arms. 



5. Conclusions 
5.1. Summciry 

A key unsolved problem in the dynamics and evolution of spiral galaxies lies in un- 
derstanding the origin of spiral-arm substructure, including gaseous spurs, giant molecular 
clouds, and complexes of giant H II regions. Various observational and theoretical arguments 
increasingly support the notion that this arm substructure may originate from gravitational 
instability of diffuse gas within spiral arms. In Paper I, we studied dynamical evolution of 
a local segment of a magnetized spiral arm in a razor-thin model of a disk. We showed that 
magneto- Jeans instability (MJl) initiated within spiral arms naturally yields gaseous spurs 
extending into the interarm region. Paper 1 further showed that these spurs fragment into 
bound condensations that could potentially evolve into arm and interarm H II regions. 

The thin-disk models of Paper 1, however, tend to overestimate midplane self-gravity, 
which favors small-scale MJl modes. The models of Paper 1 were also unable to capture 
potential dynamical consequences of the Parker instability and other three-dimensional in- 
stabihties that rely critically upon the vertical dimension. Here, we have extended Paper I to 



''Since in axisymmetry there is no stable shock for the parameters of model HU3dl, we cannot measure 
cloud spacings in terms of Asp (which is undefined). 



-27- 



consider these important effects. We consider two sets of numerical models: two-dimensional 
disks in which a thick-disk gravitational kernel (eq. [11]) approximately treats the geometric 
dilution of self-gravity, and full three-dimensional disks in which all fluid variables are al- 
lowed to vary with the vertical coordinate. As in Paper 1, both models adopt an isothermal 
equation of state. Our main objectives were to explore how the properties of spurs and 
bound condensations produced by MJI in vertically-extended disks differ from the models of 
Paper I, and to examine the effectiveness of three-dimensional dynamical instabilities other 
than MJI in forming spiral-arm substructures. The following summarizes the main results 
of the present work: 

1. In our new models with suflicicntly- magnetized conditions, spiral shocks give rise to 
gaseous spurs and bound clouds in a manner similar to the models of Paper I. This is true 
either for two-dimensional "thick disk" or three-dimensional models. The spur structures 
themselves are slightly more trailing than those in razor-thin disks. Dilution of self-gravity 
due to flnitc disk thickness is significant, causing the typical separation of spurs Xy (along 
the arm) to be about 10 times the Jeans length Aj^sp = c^/lC^sp) at the spiral arm peak, 
which is 3 - 5 times larger than the prediction of Xy/Xj^sp from thin-disk models. The 
agreement between three-dimensional and two-dimensional "thick disk" models implies that, 
in this strong- B-field case, the mechanism behind spur and clump formation is the MJI. 
Reduced gravity in thick disks also causes the bound condensations that form when spurs 
fragment to be more massive than in thin-disk models. The average mass of these clouds 
is ~ (1 — 3) X 10^ Mq corresponding to 6 - 12 times the thin-disk Jeans mass at the arm 
peak, and comparable to the thick-disk Jeans mass at the arm peak. The clump masses 
are comparable to the thin-disk Jeans mass at the mean unperturbed surface density of the 
disk. These bound clouds are magnetically supercritical with the mean mass-to-flux ratio 
of ~ (1 — 3)6"^^/^, and undergo significant loss of angular momentum (80% of the initial 
galactic value) via magnetic braking. 

2. Before gravitational instability sets in, our three-dimensional models exhibit time- 
dependent behavior in their spiral shock structure that is quite unlike the rapid approach to 
steady state that characterizes two-dimensional models. The three-dimensional distributions 
can be averaged over the y direction (parallel to the spiral arm) to yield an XZ shock profile. 
In these XZ profiles, the shock front generally moves back and forth relative to its mean 
position (in quasi-radial coordinate x). The flapping period of the XZ spiral shock in an 
unmagnetized model is ~ 0.3iorb, comparable to the epicyclic period at the arm peak. The 
amplitude of the flapping motion tends to be larger at large z in a given model, with radial 
excursions of ~ 0.06L^ at \z\/Hq = 1 in the /Sq — 1 model. This amplitude is a factor 3 or 
2 times larger than in models with /So = 10 or unmagnetized (/?o = oo) models, respectively. 
This flapping and other nonsteady motions do not seem to strongly affect development 
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of gravitational instabilities, however, probably because they depend on a mean density 
enhancement over large scales that need not be strictly coherent for growth. 

3. It is informative to consider the time average of the XZ shock profiles. The shock 
front in this mean profile is in general curved in the x-z plane. In these XZ profiles, spiral 
arm regions are generally thinner than inter arm regions; the ratio of scale heights in arm 
to interarm regions is about 0.5 for models with /3o = 10 or oo and ~ 0.8 for /3o = 1 
models. The surface density distributions of the time-averaged profiles are similar to those 
of one-dimcnsional steady spiral shocks resulting from two-dimensional thick-disk-gravity 
models averaged over y. This suggests that the thick-disk gravitational kernel of equation 
(11) provides a good approximation for the gravitational potential near the midplane of a 
three-dimensional disk. It also shows why, due to the right average conditions, gravitational 
instabilities grow, despite significant stochasticity in the fiow. 

4. We find that weakly magnetized or unmagnetized two-dimensional models are un- 
stable to the "wiggle instabihty" described by Wada & Koda (2004) (based on non-self- 
gravitating models with strong shocks). Wada & Koda (2004) advocated Kelvin-Helmholtz 
instabilities as an explanation for the wiggle instability, and argued that spiral arm spurs 
could arise from this mechanism. Indeed, we find the magnetic field topology and the gen- 
eration of vorticity near the shock front in our two-dimensional unmagnetized or weakly 
magnetized models are suggestive of Kelvin-Helmholtz instabilities in shear layers. We also 
find that mergers of vortical clumps in these two-dimensional models can eventually produce 
self-gravitating clouds with mass ~ 1 x W Mq, although the spur structures are much less 
prominent than in our two-dimensional models with stronger magnetic fields. 

Most importantly, however, we find that the vorticity-gcncrating wiggle instability is 
absent in full three-dimensional models of all magnetic field strengths. It appears that radial 
flapping motions of the XZ shock front, combined with strong vertical shear of horizontal 
velocities, quickly disrupt any coherent vortical structures that would otherwise grow. We 
therefore conclude that the wiggle instability is an artifact of two-dimensional models within 
a certain parameter range, and is unlikely to play an important role in forming spiral-arm 
substructures in real spiral galaxies. 

5. While the Parker instability has long been invoked as a primary mechanism for 
the formation of giant molecular clouds inside spiral arms, our three-dimensional models do 
not show any noticeable evidence of developing Parker modes. There is no indication of 
sinusoidal vertical velocities, or correlation of magnetic hills/valleys with over/under dense 
regions. Like the wiggle instability, the Parker instability appears to be suppressed by strong 
vertical shear of inplane velocities. 
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6. In addition to MJI, two other mechanisms capable of forming gravitationally bound 
clouds in disk galaxies are swing amplification and collapse of spiral arms parallel to the shock 
front (followed by fragmentation). The simulation outcomes of our unmagnetized, three- 
dimensional models with spiral arms show, however, that the growth of nonaxisymmetric 
disturbances via swing amplification is very small when F = 5% and Qq = 1.8. This 
is consistent with the results of Paper I, in which we showed that swing amplification in 
thin disks is efficient only if the background spiral perturbation does not exceed 1%. With 
reduced self-gravity, swing amplification in three-dimensional spiral arms becomes even more 
inefficient. Thus, the swing mechanism would appear to apply either in interarm regions or 
in galaxies without strong spiral arms. 

On the other hand, our unmagnetized model with F — 5% and Qq = 1.5 was sufficiently 
unstable that quasi-axisymmetric growth developed as soon as the external potential reached 
full strength. Fragmentation into closely-spaced bound clouds then occurred. Conceivably, 
this kind of cloud formation process could be important in galaxies that are intrinsically 
quite close to instability, and experience a tidal encounter that tips them "over the edge" . 

5.2. Discussion 

The summary above makes clear that our new three-dimensional models uphold the 
principal conclusions of Paper I - namely, that gravitational instability in the gas component 
of spiral arms (1) is able to form bound clouds with properties similar to the most massive 
GMCs, and (2) can lead to the development of interarm extensions similar to features that 
have been described in the observational literature as spiral arm spurs and feathers. 

Our models show that gaseous spurs may continue to stand out against the interarm 
background for a long time after clouds form, without being dispersed (see Figures 10 and 
11). Since interarm regions are characterized by low gas density and strong shear, they 
are normally thought of as being hostile for gas to condense gravitationally. Our model 
simulations predict, however, that even interarm regions with low mean gas density can be 
abundant with kpc- and larger-scale substructure, arranged in traiUng gaseous spurs that 
may host H II regions. The material that makes up these interarm spurs consists of parcels 
of gas that were the center of attraction during the most recent spiral arm transit. The 
concentration of interarm gas into secondary spurs may be crucial in enabling clouds, and 
hence stars and H II regions, to form at large distances from the gas-dust ridges that mark 
the primary loci of spiral arms. Because of this, azimuthal offsets between primary dust 
lanes and maxima in Ha emission may not give meaningful measures of star formation "lag" 
timescales (cf. Mouschovias, Tassis, & Kunz (2005)). 
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While the bound clouds produced in our models continue collapsing in a runaway fash- 
ion, our numerical resolution is inadequate to follow their evolution beyond the point where 
the central cloud density has grown by more than two orders of magnitude. This is because 
the local Jeans length must be resolved with at least a few spatial zones for accurate evolu- 
tion; unphysical fragmentation may otherwise occur (Truelove et al. 1997, 1998). Although 
the future of the massive, bound clouds we identify cannot be known at this point, we specu- 
late that their internal turbulent dynamics (unresolved in the present models) would lead to 
fragmentation into GMCs, some portions of which would undergo clustered star formation. 
The eventual result would be giant H II regions in both spiral arm and interarm regions. 

The masses, magnetic fluxes, and angular momenta of the clouds formed via gravita- 
tional instability in the present three-dimensional models are quite similar to the results of 
Papers II and III, in spite of differences in the specific mechanisms involved in each case. 
In Paper II, there was no spiral structure present, while in Paper III there was strong tur- 
bulence associated with the magnetorotational instability. Thus, wc conclude that it would 
be difficult to discern the detailed route to cloud formation by observing the end result; 
not enough information is retained. On the other hand, the similarity of these bulk cloud 
properties from models with rather different particulars indicates a robustness of conclusions 
for the outcome of self-gravitating cloud formation in all kinds of disk galaxies. 

Our flnding that gas in spiral arms is thinner than in interarm regions appears to be 
consistent with Nakanishi & Sofuc (2003) who reported that the scale height of the H I layer 
in our Galaxy systematically decreases with increasing midplanc density. The fact that the 
Milky Way molecular gas, which predominantly traces spiral arms, has about half the scale 
height of H I gas (Malhotra 1994) also supports our result. Stark & Lee (2005a) furthermore 
find that massive CO GMCs, which are preferentially associated with spiral arms in the 
Milky Way (Solomon, Sanders, & Rivolo 1985; Solomon & Rivolo 1989), have a smaller scale 
height than low-mass molecular clouds. We note however that more quantitative comparison 
would require consideration of ISM turbulence and other processes that may vary spatially, 
and significantly affect the gas scale height. 

In the present models, we find that 20-30% of the total gas has collected into massive 
bound clumps by the end of the simulation, in those cases where gravitational instability 
occurs. We have found similar, or slightly smaller fractions in the uniform-shear (go = 1) 
three-dimensional models of Paper II and 111 (the range is ~ 15 — 20%). This fraction can be 
thought of as the efficiency, per GMC formation timescale, of converting diffuse gas to gas 
capable of forming stars. In the present models, the timescale for gravitational instability to 
run away is 2-6 orbits, for a range of values of Qo, F, and /?o- Combining these raw numbers 
would imply an effective GMC formation timescale from diffuse gas of ~ lO^orb- However, 
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we believe that in fact this somewhat overestimates the true timescale for the average parcel 
of diffuse gas to be incorporated in a bound cloud. While the conversion efficiency per 
GMC formation timescale may be realistic, the time to reach runaway is probably longer 
than it would be in a real system, because our initial perturbation amplitudes arc likely low 
compared to realistic values. An important direction for future study will be to measure 
cloud formation rates in simulations where cloud destruction is also incorporated; we plan 
to pursue studies of this kind. Models of this kind will provide more realistic background 
conditions from which gravitational instabilities can develop. 

While the real ISM has a multi-phase structure due to the particular heating and cool- 
ing processes that are relevant, we have adopted a much simpler isothermal approximation 
throughout this paper. The inclusion of gaseous cooling that leads to dense cloud formation 
at small spatial scales (e.g., Koyama & Inutsuka 2002; Piontek & Ostriker 2004, 2005; Audit 
& Hennebelle 2005; Heitsch et al. 2005; Vazquez- Semadeni et al. 2005b) might indirectly 
affect formation timescales, mean separations, and masses of the GMCs and other structures 
that result from self-gravitating instabilities in spiral arms. It is often assumed that the to- 
tal velocity dispersion in the warm/cold medium can be treated as an effective sound speed; 
whether this treatment is adequate will be addressed directly in future work. Our three- 
dimensional isothermal simulations show that the wiggle instability is stabilized by strong 
vertical shear and non-steady flows associated with spiral shocks, but the situation could 
conceivably be different if thermal instability and multi-phase gas structure are considered. 
This outstanding issue will direct our future research as well. 

Finally, we remark that, while spiral structure clearly can help gravitational instability 
and prompt star formation, it may also contribute to limiting these processes. As we have 
discussed, spiral shocks (especially when strongly magnetized), lead to time dependent mo- 
tions that may be important in feeding turbulence at scales <^ H. The present models are 
isothermal at T ~ 8, OOOK, so the sound speed 7 km s^^ exceeds turbulent velocity ampli- 
tudes. However, in the real ISM there is significant cold gas for which the random turbulent 
motions far exceed the sound speed ~ 1 km s~^. The turbulence in the cold ISM component 
may be key to regulating self-gravitating instabilities. In a future publication, we will discuss 
results of an investigation focusing on the details of turbulent driving by spiral shocks. 
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Table 1. Parameters of Two-Dimensional Models with Thick-Disk Gravity 



ModeP 


Qo 




F 


Resolution 




Qsp 






Ky 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


ME2dl 


1.2 


1 


10% 


256 X 512 


4.20 


0.59 


0.15 


0.50 


0.11 


ME2d2 


1.2 


1 


10% 


128 X 256 


4.14 


0.59 


0.15 


0.50 


0.11 


ME2d3 


1.2 


1 


10% 


128 X 128 


4.14 


0.59 


0.15 


0.67 


0.08 


MS2dl 


1.5 


10 


5% 


256 X 512 


5.96 


0.61 


0.062 


0.22 


0.22 


MS2d2 


1.5 


10 


5% 


128 X 256 


5.53 


0.64 


0.069 


0.22 


0.23 


MS2d3 


1.5 


10 


5% 


128 X 128 


5.53 


0.64 


0.069 


0.33 


0.16 


H2dl 


1.8 


00 


5% 


256 X 512 


6.75 


0.69 


0.044 


0.16 


0.32 


H2d2 


1.8 


00 


5% 


128 X 256 


6.33 


0.72 


0.049 


0.20 


0.27 


H2d3 


1.8 


oo 


5% 


128 X 128 


6.33 


0.72 


0.049 


0.31 


0.18 



''The prefixes ME, MS, and H stand for the magnetized models with equipartition (/3o = 1) and sub-equipartition 
(/3o = 10) field strengths, and hydrodynamic models. 
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Table 2. Parameters of Three- Dimensional Simulations 



Model 


Qo 




F 


Resolution 






-ffarm/ -^0 




Ky 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


ME3dl 


1.2 


1 


10% 


128 X 256 X 128 


3.32 


0.66 


0.93 


0.50 


0.14 


ME3d2 


1.2 


1 


10% 


128 X 128 X 128 


3.30 


0.66 


0.93 


0.67 


0.10 


MS3dl 


1.5 


10 


5% 


128 X 256 X 128 


5.31 


0.65 


0.71 


0.50 


0.11 


MS3d2 


1.5 


10 


5% 


128 X 128 X 128 


5.30 


0.65 


0.72 


0.67 


0.08 


H3dl 


1.8 


00 


5% 


128 X 256 X 128 


5.06 


0.80 


0.64 






H3d2 


1.8 


oo 


5% 


128 X 128 X 128 


5.02 


0.80 


0.66 






HU3dl^ 


1.5 


oo 


5% 


128 X 256 X 128 








0.15 


1.92*^ 



'^The prefix HU stands for hydrodynamic model with a strongly unstable arm. 
'^Relative to the initial background disk without spiral potential perturbations. 
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Fig. 1. — Time evolution of maximum surface density for two-dimensional high- resolution 
models. 
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(b) t = 2.5t,,. 




(c) t = 3.2t„ 
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0.5 -0.5 
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0.5 



Fig. 2. — Snapshots of model ME2dl at t/torb = 2.5 and t/torb = 3.2. For fiducial parameters, 
the simulation box has a size of = Ly/2 = 3.14 kpc. (a) Surface density and velocity 
fields seen in the comoving frame with the spiral pattern, (b) Surface density and a few 
selected streamlines {red) in the frame comoving with the spurs, together with wavefronts 
{black) defined by equation (12) of Paper I with Kx{0) = 0.45. The ?/-boundaries are shifted 
so as to make the left wall coincide with the shock location, (c) Magnetic field lines {red) 
are overlaid on surface density. Colorbars label logS/So. 
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Fig. 3. — Prediction of spiral arm MJI development from linear theory, with thick-disk 
gravity, for parameters of model ME2dl. The gaseous scale height Hq — 180 pc is taken 
as fixed. Solid contours, spaced at 0.4, 0.6, 0.8, • • • , 1.6, 1.8 from outside to inside, show final 
amplification magnitude F = log I^Emax/^^l as a function of Ky and Kx{0), the normalized 
wavenumbers at the shock location parallel and perpendicular to the spiral arm, respectively. 
Dotted contours plot the corresponding growth times tgrow/^orb = 0.2,0.4,0.6,0.8,1.0 from 
right to left. The most unstable mode evident in the simulation results from model ME2dl 
is indicated by the box near Ky = 0.11 and Kx{0) = 0.45. 
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Fig. 4. — Dependence on the gaseous scale height Hq of (sohd hnes; left y-axis) the wavenum- 
bers i^a:,max(0) and i^jy,max and (dotted line; right y-axis) the amplification magnitude Fmax 
of the dominant MJI mode in spiral shocks predicted from linear theory. The parameters 
for model ME2dl are used. 
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Fig. 5. — Surface density maps (log(S/So) in gray scale) of model MS2dl at (left) t/torb = 1-6 
and (right) t/to^\^ = 2.7. For fiducial parameters, the simulation box has a size of = Ly/2 = 
3.14 kpc. In both panels, the initial shock location is indicated by the dotted line, while 
the dashed line marks the position where the sense of shear in the initial background flow 
changes from reversed to normal. The rectangular section in the left panel is enlarged in 
Figure 6 to show its temporal changes. 
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Fig. 6. — Temporal changes of the rectangular section shown in the left panel of Figure 5 
at (a) t/torb = 1-1, {b) t/torh = 1-3, (c) t/torb = 1-6, and {d) t/torh — 1-9. In each panel, 
magnetic field lines are overlaid on surface density (log(E/Eo) in gray scale). 
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Fig. 7. — Maximum gas density in units of the initial midplane density po(0) for high- 
resolution three-dimensional models. The imposed spiral potential reaches full strength at 
t/toTh = 1-5. While the density in unmagnetized model H3dl fluctuates around a mean 
value and does not become sufficiently nonlinear for gravitational instability, magnetized 
models ME3dl and MS3dl, which are subject to MJI, evolve into a highly nonhnear state, 
forming spurs and gravitationally bound condensations. Unhke models MESdl, MS3dl, 
and HSdl which permit stable shock equihbria, model HU3dl produces a spiral shock that 
grows exponentially via quasi-axisymmetric Jeans instability and subsequently experiences 
non-axisymmetric fragmentation to form bound condensations at the shock front. 
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Fig. 8. — Mean spiral shock solutions in the radial-vertical (XZ) plane of model MS3dl, 
based on spatial average along azimuth {y) and time average over t/torb = 1.7 — 2.7. (a): 
Gas streamhnes (thin sohd hues) are overlaid on the gas density (log(p/po) in gray scale) . 
The heavy vertical hue marks the averaged shock front, (b): Average surface density profile 
(solid line) based on vertical integration of density in (a) is compared to its two-dimensional 
thick-disk counterpart from model MS2d2 (dotted line), (c): The profile of mean gaseous 
scale height for models MS3dl (solid line) and MS2d2 (dotted line). The gaseous scale height 
is generally smaller inside spiral arms than in inter arm regions. 
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Fig. 9. — (left) Vertical distributions and (right) gradients of spiral shock velocities. Mea- 
surements are slightly upstream from the midplane shock at x/L^ — —0.05 (dotted), slightly 
downstream from the shock at x/L^ — 0.05 (solid), and in the interarm region aXx/L^ = 0.22 
(dashed). Profiles are based on the averaged XZ shock configuration of model MS3dl. Shear 
in the horizontal velocities is very strong for 0.2 ^ \z/H\ ^ 2. 
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Fig. 10. — Snapshots of vertically-integrated surface density (log(S/Eo) in gray scale) of 
(3o = 10 model MS3dl at (a) t/torb = 5.6, (6) t/torb = 6.0, and (c) t/torb = 6.3. For fiducial 
parameters, the simulation box in the x-y plane has a size of — Ly/2 — 3.14 kpc. The 
rectangle in (a) indicates the sector viewed as a three-dimensional visualization in Fig. 12. 
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Fig. 11. — Snapshots of surface density (log(S/So) in gray scale) of /3o = 1 model ME3dl at 
(a) t/torb = 3.0, (6) t/torb = 3.1, and (c) t/torb = 3.2. For fiducial parameters, the simulation 
box has a size of — Ly/2 — 3.14 kpc. Due to the strong radial flapping motions of gas 
near the spiral shock, spurs in model MESdl are less conspicuous than in model MS3dl. 
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Fig. 12. — (a) Volumetric rendering (seen from above the galactic plane) of an isodensity 
surface at p = 2.5po(0), together with selected magnetic field hues, from a portion of model 
MS3dl at t/torb — 5.6. The region fits in a projected rectangle marked in Fig. 10, with 
\z\/ H <1.25. The magnetic field lines in blue run from {y, z) = {0.35Ly, 0) at the lower edge 
of the box toward the upper edge, while yellow denotes field lines originating at {y, z) = 
{0.35Ly,0.56H). The density ((log(p/po) in color scale) and velocity vectors at z = are 
shown at the bottom of the box, while those at y = 0A5Ly are displayed on the side wall 
(top of page), {b) The density (log(p/po) in color scale), velocity vectors, and magnetic field 
lines from the left surface {x/Lx — 0.5) of the box shown in (a). In both panels, the velocity 
vectors are relative to the center of the box, where the density attains its maximum. 
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Fig. 13. — YZ slice in the upper half plane at x/ = 0.05 of model MS3dl at t/torb = 4.62 
when spurs are about to emerge. The density field (log(p/po) in color scale) , perturbed 
velocity vectors, and magnetic field lines are drawn. Strong vertical shear of the azimuthal 
velocity is apparent. The characteristics of the Parker instability such as the correlation 
between density and magnetic fields as well as the sinusoidal variations in the vertical velocity 
are not evident. 
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